%--------------------------------------------------------------------------
%Smoothing Quadratic Regularization Methods for Non-Lipschitzian
%Optimization

%Solve
%          min   f(x):= H(x)+ \sum_{i=1}^n \phi(|x_i|^p), x \in R^n
%                                                                     
%Please refer to Test_SQR to get information on how to use this algorithm
%
%Reference: "Worst-case Complexity of Smoothing Quadratic Regularization
%            Methods for Non-Lipschitzian Optimization"
%            W. Bian and X. Chen, SIAM J. Optim. 23 (2013), 1718-1741
%                                                                                                             
%--------------------------------------------------------------------------

%--------------------------------------------------------------------------
%Explanation of Parameters:
%    Input Parameters:
%    x0: the initial point of iteration, it's a column vector
%    mu0: the initial value of smoothing parameter
%    beta0: a constant that is not smaller than 1, used to determine the
%           value of regularization weight
%    beta1: a constant that may be the upper bound of the Hessian of H(x)
%           if H(x) is second differentiable
%    sigma: 0<sigma<1, a constant, used to update the smoothing parameter
%    sigma1: 0<sigma1<1, a constant, used to update the regularization
%            weight
%    sigma2:0<sigma2<1, a constant, used to judge whether the iteration is
%           very successful or succussful
%    eta: eta>1, a constant, used to update the regularization
%    alpha: the common bound of $\phi\prime$,$|\xi|$,$|\xi|t, \forall \xi \in
%           \partial (\phi\prime(t)), \forall t \in (0,\infty)$
%    p: 0<p<=1, a constant appears in the objective function
%    eps: the $\epsilon$ of the definition of $\epsilon$ scaled stationary
%         point

%    Output parameters:
%    We use a struct named Out to store the output paremeters.
%    Out.x: the iteration sequence of x
%    Out.z: the convergent sequence z
%    Out.sgr: the gradient sequence of smoothing function respecting to the 
%             convergent sequence z
%    Out.xsgr: the product sequence of z and its gradient of smoothing function 
%--------------------------------------------------------------------------


function Out=SQR_NonLip(x0,mu0,beta0,beta1,sigma,sigma1,sigma2,eta,alpha,p,eps)

n=length(x0);
Gtol=eps;
C=8*alpha*p;


% Initialization
xk=x0;
zk=xk;
mu=mu0;
beta=max(beta1,beta0);

iter=1;
Rec_iter(iter)=iter;
x(:,iter)=xk;
z(:,iter)=zk;

xsg=sgrad(xk,mu,p);
zsg=sgrad(zk,mu,p);
GTsg=sgrad(zk,Gtol,p); %Used to judge if the initial point is the minimal point

xsgrad(:,iter)=xsg.gr;
xxsgrad(:,iter)=xsg.xgr;
zsgrad(:,iter)=zsg.gr;
zzsgrad(:,iter)=zsg.xgr;


% Main Loop
fprintf('\n Iteration:    ');

while (zsg.xgr>Gtol || mu>Gtol) && GTsg.xgr>Gtol
    fprintf('\b\b\b\b%4d',iter);
    %compute the step
    step_length=cal_step(xk,mu,beta,p);
    %new point calculation
    yk=xk+step_length;
    
    %compute r_k
    if norm(yk-xk)==0
        rk=1;
    else
        rk=(H(yk)-H(xk)-grad_H(xk)'*(yk-xk))/(0.5*beta*norm(yk-xk)^2);
    end
    
    %Update the regularization weight
    if rk>1
        beta=eta*beta;        
    elseif rk<sigma2
        beta=max(beta0,sigma1*beta);        
    end
    
    %Update the smoothing parameter and construct convergent point
    if rk<=1
        xt=xk;
        xk=yk;
        if sfunc(xk,mu,p)-sfunc(xt,mu,p)>-0.5*C*mu^p
            mu=sigma*mu;
            zk=xt;               
        end
    end
    
   
    xsg=sgrad(xk,mu,p);
    zsg=sgrad(zk,mu,p);
    
    iter=iter+1;
    Rec_iter(iter)=iter;
    % Record data
    x(:,iter)=xk;
    z(:,iter)=zk;
    
    xsgrad(:,iter)=xsg.gr;
    xxsgrad(:,iter)=xsg.xgr;
    zsgrad(:,iter)=zsg.gr;
    zzsgrad(:,iter)=zsg.xgr;
     
end

Out.x=x;
Out.z=z;
Out.sgr=zsgrad;
Out.xsgr=zzsgrad;
Out.iter=Rec_iter;




%--------------------------------------------------------------------------
%SUBFUNCTION FOR CALCULATING smoothing function, gradient of smoothing
%fuction and step
%--------------------------------------------------------------------------

%smoothing function of the objective function
    function sf=sfunc(x,mu,p)
        
        sf=H(x)+sum(phi(theta(x,mu).^p));       
    end

%gradient of smoothing function
    function sg=sgrad(x,mu,p)
        
        idx1=find(abs(x)>mu);
        idx2=find(abs(x)<=mu);
        
        dthetap(idx1,1)=p*abs(x(idx1)).^(p-1).*sign(x(idx1));  
        x_dthetap(idx1,1)=p*abs(x(idx1)).^p;
        dthetap(idx2,1)=p*(x(idx2).^2/(2*mu)+mu/2).^(p-1).*x(idx2)/mu;
        x_dthetap(idx2,1)=p*(x(idx2).^2/(2*mu)+mu/2).^(p-1).*x(idx2).^2/mu;
        t=theta(x,mu).^p;
        dphip=Dphi(t).*dthetap;
        xdphip=Dphi(t).*x_dthetap;
        
        sg.gr=grad_H(x)+dphip;
        sg.xgr=max(x.*grad_H(x)+xdphip);
        
    end


%calculate step of iteration
    function h=cal_step(x,mu,beta,p)
        sgra=sgrad(x,mu,p);
        sg=sgra.gr;
        abs_sg=abs(sgra.gr);
        
        temp1=mu./abs_sg;
        temp2=1/(beta+C*mu^(p-2));
        Rec_gamma1=min(temp1,temp2);
        
        temp3=mu^(p/2)*abs(x./2).^(1-p/2)./abs_sg;
        temp4=1./(beta+C*abs(x./2).^(p-2)); 
        Rec_gamma2=min(temp3,temp4);
        
        Index=zeros(n,1);
        Ones=ones(n,1);
        idx=find(abs(x)>2*mu);
        Index(idx,1)=Ones(idx,1);
        Rec_gamma=(1.-Index).*Rec_gamma1+ Index.*Rec_gamma2;
        
        h=-sg.*Rec_gamma;
    end

%smoothing function of the absolut value,|*|
    function th=theta(x,mu)
        idx1=find(abs(x)>mu);
        idx2=find(abs(x)<=mu);
        th(idx1,1)=abs(x(idx1));
        th(idx2,1)=x(idx2).^2/(2*mu)+mu/2;
    end
        
          
end








    
  













